import numpy as np
from scipy.integrate import quad, simps, cumtrapz
from rcsj_model_analysis import inference222 as infer
from rcsj_model_analysis import model
from unittest import TestCase
from scipy import constants
h_bar = constants.hbar
k_B = constants.Boltzmann
elementary_charge = constants.elementary_charge


class TestInference(TestCase):

    def test_hazard_func(self):
        "tests the hazard function"
        I = np.linspace(0, 0.08, 100)
        critical_current = 0.08
        f_j0 = 10**9
        Temp = 30
        
        current = I
        eta = current/critical_current
        freq_j = f_j0*(1-eta**2)**(1/4)
        E_J = h_bar*critical_current/(2*elementary_charge)
        Delta_U = 2*E_J*(np.sqrt(1-eta**2)-eta*np.arccos(eta))
        test_array1 = freq_j*np.exp(-1*Delta_U/(Temp*k_B))
        
        array1 = infer.hazard_func(I, critical_current, f_j0, Temp)
        
        self.assertEqual(array1.shape,test_array1.shape,'Shape does not match')
    
    def test_loglike(self):
        "tests the loglike function"
        I = np.linspace(0, 0.08, 100)
        critical_current = 0.08
        f_j0 = 10**9
        Temp = 30
        numpoints=100
        theta = [critical_current,f_j0,Temp]
        s = np.linspace(0, I, numpoints)
        integrals = simps(infer.hazard_func(critical_current, f_j0, Temp, s), s, axis=0)

        test_array2 = np.sum(np.log(infer.hazard_func(critical_current, f_j0, Temp, I)) - integrals)

        array2 = infer.loglike(theta,I,numpoints=100)

        self.assertEqual(array2.shape,test_array2.shape,'Shape does not match')


        
if __name__ == '__main__':
    unittest.main()
